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INTRODUCTION 


Formulae  for  the  calculation  of  the  estimated  power  spectrum 
and  its  variance  from  the  true  value  for  a stationary,  random  process 
given  a finite  record,  have  been  given  by  Tukey.  These  methods 

have  been  applied  by  the  author  to  determine  the  power  spectra  of 
temperature  records.  Since  the  IBM  650  computation  programs  (written 
by  A.  Anastasia)-  for  the  evaluation  of  these  quantities  may  be  gener- 
ally useful,  they  will  be  described  in  detail,  The  mathematical 
background  will  be  developed,-) with  the  notation  and  terminology  of 
Reference  3*?  and  some  of  the  different  notations  encountered  in  the  ex- 
tensive literature  in  this  subject  will  be  pointed  out. 


THEORY 


We  consider  a function 


GT(t)  = G(t) 
= 0 


-T  < t < T 
outside. 


(1) 


The  Fourier  transform  of  G^,(t)  is 


uu  , 

S(f)  = j*  GT(t)e"lf,'tdt 

«J 


0^ 


= 2tt  f 


(2) 


and 


GT(t)  = j S(  f)ei°'tdf 


(3) 


I 
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Parseval's  Theorem  states  that 


OU  "P  ou 

J |GT(t)|2dt  = J | G( t ) | 2dt  = J |S(f)|2df 


(4) 


so  that  the  time  average  of  G (t)  is 


G2(t)  s lim  ~ [ | G( t ) J 2dt 
_JT 
OO 

= lim-^  r |S(f)|2df 

- OO 

OO 

= J P(f)df 


(3) 


where  P(f)  = lim  |S(f)|2 


T-*-  > 


(6) 


P(f)df  represents  the  contribution  to  the  time  average  of  G^it)  from 
frequencies  between  f and  f + df  and  is  called  the  power  spectrum  or 
spectral  density  or  power  spectral  density-  Note  that  P(f)  so  defined 
is  a two-sided  function  of  frequency  whereas  common  practice  is  to  use  the 
one-sided  power  spectrum,  Q(f)  . That  is, 

Q(f)  = 2P(f)  (7) 

which  is  the  power  density  for  positive  frequencies.  The  utility  of  using 
the  two-sided  function  lies  in  the  simplification  of  certain  integrals. 

We  define  the  autocovariance  function  (also  called  the  auto- 
correlation function  in  the  literature)  to  be 


l T 

C(x)  = lim  ~ P G(t)G(t+T)dt 

T->«  i 
-T 


(8) 
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where  T is  the  time  lag.  The  variance  is  C(0)  . The  Fourier  transform 
of  C(x)  is  the  power  spectrum. 


P(f)  = C(T)e'i0TCdx 


which  follows  if  one  substitutes  the  Fourier  transform  of  G(t)  in  the 
integral  (8)  and  conversely, 


C(t)=J  P(f)eia'n:df  , 


Since  C(x)  is  an  even  function  of  x . the  power  spectrum  is  also 


P(f)  = 2 J C(x)  cos  anrdr 


The  one-sided  power  spectrum  is 


Q(f)  = 4 j^CXx)  cos  amdr 


which  is  a form  quite  often  encountered  in  the  literature. 

In  practice,  our  records  are  rather  short  and  we  would  like  to 
do  the  processing  digitally  thereby  avoiding  the  necessity  of  constructing 
apparatus  for  the  evaluation  of  the  integrals  (8)  and  (9)  and  the  larger 
errors  inherent  in  analog  processing. 

For  numerical  processing,  the  data  should  be  read  at  equal  time 
intervals.  Let  the  N+  1 numbers 

V V ¥2 rN  (13) 

be  the  time  series,  taken  at  time  intervals  At  , with  mean 

hi.1,  • (14’ 
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X = Y - Y 

q q 


be  the  fluctuation  from  the  mean.  We  evaluate  the  integral  (8)  as  the 
mean  lagged  product,  with  the  maximum  lag  m 


n - r 

; = Y H ^ r = 0,1,2,  ...  m . 

r N - r qq+r 


Then,  the  Fourier  transform  (corresponding  to  Eq.  (11)  is  the  raw  spectral 
density  estimate 


V = At  C + 2 Y C cos  •2II4  C cos  m 
r o q mm 


where  is  the  estimate  of  the  two-sided  spectral  density  centered 

about  the  frequency 


f = (r/2mAt)  r = 0,1,2,  ...  m 


Factoring  the  variance,  Cq  , we  find 


where 


V = C AtV  ' 
r o r 


m - 1 

V ' - 1 + 2 Y C 1 cos  -asr+c  ' cos  rrr 
r q mm 

q=  i 


c 1 = (c  /c  ) . 
q q 0 


The  one-sided  power  spectrum  estimate  Wf  is 


W = 2V 
r r 
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We  note  that  the  sum 


t »r  Af  = 2Ct  Z vr 

r=  0 r 


>-1, 


(22) 


= 2(2mAt)  (C  Atm) 
o 


= C 


since  Af 


4 


(2mAt)  1 for  r / 0,m 


= ^(4mAt)  * for  r = 0,m 


for  the  contribution  from  the  cosine  terms  in  the  summation  over  r is 
zero. 

An  alternate  formulation  is  to  refer  the  power  spectrum  estimate 
to  unit  frequency  bands,  i.e.,  Af=  1 . That  is,  let 

(23) 


then. 


L = V (l/2mAt) = (C  /2m)V  ' , 

r r or 


£ 2Lr  = 2(l/2mAt)£Vr  = 


as  before. 

The  are  the  Fourier  transform  of  the  mean  lagged  products 
so  that  formally,  the  problem  is  solved.  However,  the  are  subject  to 

the  errors  arising  from  the  mathematical  approximation  in  their  evaluation 

as  well  as  to  the  complex  propagation  of  errors  of  the  original  data.  To 

make  a direct  computation  of  these  errors  involves  considerable  effort;  at 

least  as  great  as  for  the  computation  of  the  V . 

We  consider  instead  the  variability  to  be  expected  in  the  power 
spectrum  of  a stationary,  random  Gaussian  time  series.  For  example,  con- 
sider a large  number  k of  oscillators,  each  with  a different  frequency, 


with  random  phases  and  with  amplitudes  in  a Maxwellian  distribution.  The 

2 2 

energy  radiated  per  oscillator  is  (a  +b  ) where  the  amplitude  of  the  rad 
ation  is  given  by  (a  cos  c»t+b  sin  oit).  The  sum  of  the  squares  of  the  2k 
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numbers,  I]  |a^ +b^J  is  the  rate  at  which  the  ensemble  energy  is  radi- 
ated. This  sum  has  a Chi-squared  distribution  with  2k  degrees  of  freedom. 
Hence,  if  the  power  spectrum  density  in  any  given  frequency  band  depends 
on  v degrees  of  freedom,  we  can  use  tables  of  Chi-squared  to  determine 
the  variability  of  the  calculated  power  spectrum  from  the  true  value.  While 
not  strictly  correct,  this  result  gives  a good  estimate  of  the  variability 
of  the  power  spectrum  The  number  of  degrees  of  freedom  associated  with  the 

Vr  ,s  v-T'2  ' 

TABLE  I 

4 

Reliability  of  power  spectrum  estimates 


Degree  of  freedom 

2.5% 

Possible 

5% 

error 

95% 

97.5% 

1 

1000 

250 

0.26 

0,2 

2 

40 

20 

0.33 

0.21 

5 

6 

4.37 

0.46 

0.39 

8 

3.8 

2.8 

0.51 

0.46 

10 

3.1 

2.6 

0.55 

0.49 

15 

2.4 

2.1 

0.60 

0.55 

20 

2.1 

1.8 

0.63 

0.59 

50 

1-55 

1.45 

0.74 

0.69 

100 

1.35 

1.29 

0.79 

0.78 

The  true  value  of  the  spectrum  level  will  lie  between 
2.6  and  0.55  of  the  calculated  value  of  90  percent 
of  the  time  if  there  are  10  degrees  of  freedom  for 
the  estimate. 
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We  can  determine  the  effective  width,  or  window,  in  the  frequency 
domain  of  the  raw  spectral  density  estimates  by  computing  the  response  to 


a sine  wave. 


Let  the  input  signal  be 
X = a cos 


|27rf0t+  6 


which,  sampled  at  time  intervals  At  . becomes 


= a cos  -^P+gJ  q = 0,1,2,  ... 


t — qAt 

f =-*- 
o PAt 


ii  - r 

C<r,=^,S 

t0cos  [p'c^r26] 


The  summation  becomes 


sin  ^ (N  - r+  1) 


cos  +20 1 < 1 


so  that  we  can  ignore  it  compared  to  the  first  term.  The  raw  spectral  density 
estimate  (Eq.  17)  is 

n o ~ 

V = 2At  y C cos  q^-Atfc  +C  cos  rrr ] . 

* ntn  Q m Lorn  J 


The  first  term  becomes 


/a2\  ™ f Trqu+  irqu 

a,(t  50rs~+cos 


Trqu_1 

m J 
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where 


I 


u+  = r+  2m 
u = r - 2m 


which  reduces  to 


/TTU+,  77-U  ,,  /TTU  , TTU  , 

cos>~lsin  -rl^ml  cost— )sin  -gi^m) 


For  m large,  1-*-—  =»  1 , so 
m 


/ 2\  sin  mi  sin  w 

Vr=%) ^ + S 

2sI"  2iT  2sin  2m 


-At[VC„ 


cos  rir 


In  the  neighborhood  of  small  u , that  is,  near  the  frequency  of  the  sine 


wave, 


v ^ Jinja.  „ Ijinja 
r 2sinfJ  1 m 


which  is  1 for  m«=0  , and  has  a full-width  at  half-Hnaximum  of 

Ar  = 1*5  * M 


and  has  many  minor  lobes,  the  largest  occurring  at 


3 i 2m 

r = %*F  ' 

of  amplitude  - .2  j . This  is  20  percent  of  the  main  peak  and  is 

much  too  large. 

By  multiplying  this  filter  function  by  a factor 


■i  l+  cos  ^ 
c m 
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we  can  reduce  the  side  lobe  amplitude  to  approximately  2 percent  of  the 
main  peak.  The  response  to  a sine  wave  is  now  given  by 


.TTU, 


TTU, 


n 


2m 


r 

cos 

flv»J 

sin 

if  !•♦♦•) 

(■41. 

2sin 

k+1)l 

-1)]  sin  [f  (u+ 

->] 

1 +in  |] 

^‘"[Sk  - x)] 

+ ...  similar  terms  for 


.1 


which  may  also  be  written  as 


u 4v  ,+^V+^V  . 

r 4 r-1  2 r 4 r-1 


(24) 


These  are  called  the  refi ned  spectral  density  estimates  where  is  the 
average  power  density.  One  can  compute  that  this  average  is  taken  over  a 
filter  width  of  centered  at  the  frequency  fy  . 


ERRORS 

The  finite  approximations  involved  in  this  computation  introduce 
a number  of  sources  of  errors  if  proper  precautions  are  not  taken.  Two  of 
the  most  serious  will  be  mentioned  here. 

(a)  Aliasing  - If  the  sampling  interval  is  At  , then  the 
highest  frequency  at  which  spectral  density  estimates  can  be  made  is 
f^  = 2At  * **  appreciable  power  is  present  at  frequencies  greater  than 

f^  , they  will  appear  at  the  frequencies 

f . 2fN  * f . 4fN  * f , . . . 
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form  will 


For  example,  Fig  1 shows  how  the  aliased  spectrum  of  a f 
appear  Hence,  in  the  design  of  the  experiment,  one  must  insure  that  the 
upper  frequency  cutoff  of  the  measuring  instruments  is  less  than  f^  . 
Alternate] v,  low-pass  filtering  must  be  resorted  to  after  the  data  is  taken 
and  before  it  is  spectrum  analyzed. 

(b)  Prewhiteni ng  - If  the  true  spectrum  has  large  peaks,  then 
there  will  be  diffusion  of  power  to  adjacent  frequency  bands  through  the 
side  lobes  of  the  filter  function.  Therefore,  one  must  prewniten,  i.e.,  make 
the  spectrum  more  nearly  flat  by  suitable  processing  of  the  time  series  in 
the  time  domain,  and  correct  afterwards  in  the  frequency  domain  for  the  ef- 
fects of  prewhitening. 

An  example  of  such  prewhitening  by  either  high-pass  or  low-pass 
filtering  is  the  processing  in  the  time  domain  of  the  original  data  in  the 
form 

x. ' = x. + bx.  . -1  < b < 1 . (25) 

i i i-I 

The  effect  in  the  frequency  domain  is  to  multiply  by  a filter  function  of 
the  form 

Y2(w)  = 1+  b2  +2b  cos  w (26) 

so  that  the  calculated  power  spectrum  is 

P(calculated)  = Y2(w)P(true)  „ (27) 

In  one  of  the  programs  available,  we  have  set  b = - 0.6  The  effect  of 
this  is  to  reduce  the  low-frequency  power  by  a factor  of  0.16  relative  to 
midband,  and  to  increase  the  high-frequency  power  by  1.96  . 

Other  values  of  b can  be  used  provided  they  are  inserted  into 
the  proper  locations  of  the  program. 
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Fig.  1 


FREQ 

Fig.  2 


J 


'!!  - — 


A.  Noise 


RESULTS 


The  spectral  analysis  has  been  applied  to  a sequence  of  random 
numbers.  Figure  2 gives  the  power  spectrum  which  is  flat  with  fluctu- 
ations in  the  order  of  the  90  percent  confidence  limits  (Table  I).  Hence, 
fluctuations  of  this  order  in  experimental  data  may  also  be  due  to  random 
effects, 

B.  Signal 

The  power  spectrum  of  a CW  signal  of  the  form 
y = a(l  - cos  27rq/4) 

has  been  calculated.  The  solution  for  the  continuous  case, 

y = a cos  2?rf  t , 
o 

may  be  easily  calculated.  The  autocovariance  function  is 

C(T)  = (a2/2)  cos  27rf0T  • 

The  two-sided  power  spectrum  is 

PCD  =^[&|f+fj+6|f-f0)]  . 

and  the  one-sided  power  spectrum  is 

«<»  =CoMf-'o)  • 

The  results  of  the  Tukey  analysis  give  for  the  case  fQ  = (l/4At)  and 
m — 26  , 


2V  = 0 for  r / 13 
r 

= 2mC  At  for  r = 13  . 
o 
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raiiJirl, 


Then, 


2U  - 0 for  r A 12,  13,  14 
r 

= jmC  At/2]  for  r = 12,  14 
= mCoAt  for  r = 13  . 

The  integral  of  the  power  spectrum  over  the  frequency  domain  is  Cq  as 
it  should  be 
C Signal  plus  noise. 

A small  amplitude  CW  signal  is  detectable  in  considerable  noise 
as  the  following  analysis  shows.  Basically,  this  is  because  the  power  of 
the  CW  signal  is  concentrated  at  a single  frequency  or  in  a band  of  fre- 
quencies whereas  the  power  in  the  noise  is  spread  over  the  entire  frequency 
domain  as  was  shown  in  A. 

Let  e(t)  = noise  signal  and  a(l-  cos  cot)  = signal.  The  sum  of 
the  two  y(t)  = e(t)  + a(l-cos  is  the  function  to  be  analyzed.  Let 

the  mean  value  of  y(t)  be 

__  v 

y = e + a 
so 

x(t)  = y(t)  - y 

= e'(t)  - a cos  cnQt 

where 

e; (t)  = e(t) - e 

Then 


ou 

C(x)  = lim  -5=  [ [e ' (t)e ' (t+ x)+  cos  t cos  con(t+T) 
T-»oo  £l  J L 0 0 


-ae'(t+  t)  cos  (o  t - ae'(t)  cos  0)  (t  + x)l  dt  (C.l) 
0 0 J 


. ..  - . 
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The  first  term  in  the  integral  is  the  autocovariance  function  of  the  noise 


cn(t)  = lim  -jyj:  J e'(t)e'(t+  x)dt  . 


Similarly,  the  second  integral  is  that  of  the  signal 

1 rT  2 

as(x)  = lim  J a cos  coot  cos  <oQ(t+  x)dt 


The  last  two  terms  are  the  cross  correlation  of  the  signal  and 

noise 

<Wt) 

which,  in  the  lim  acN^J-^O  since  the  signal  and  noise  are  uncorrelated. 

T-+.OQ  biN 


We  can  rewrite  (C„l)  as 


C(x)  = aN(x)+  ffg(T)+  Csn(t) 
~ aN(Tr)+  cs(t)  . 


(C.2) 


The  Fourier  transform  of  C(t)  is 

P(f)  = PN(f)+  ps(f)+  pSN(f) 

~PN(f)+  Ps(f) 

We  can  immediately  apply  (C.3)  to  the  digital  computation. 

From  A , the  one-sided  power  spectrum  of  the  noise  is 


(C.3) 


N 

Wr  = ^ = 2 ‘NAt 


where  the  letter  N designates  the  noise  spectrum. 
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We  see  that  a CW  signal  of  small  variance  relative  to  the  noise  can  be 
detected  provided  m is  large  enough.  For  a CW  signal,  of  amplitude  a , 


°S 


so  the  minimum  detectable  signal  is 

3 = aN 

and  if  = e,  and  p = 2 , 


(C.5) 


PROGRAM  DETAILS 

Two  programs  are  currently  available  for  the  computation  of  power 
spectra.  The  first  computes  the  power  spectra  of  short  time  series  with- 
out preprocessing  of  the  data.  The  second  allows  for  the  processing  of  the 
data  by  prewhitening,  can  handle  much  longer  time  series,  but  requires  three 
steps  to  complete.  There  is  also  available  a program  for  directly  computing 
the  Fourier  series  expansion  of  a time  series  but  its  use  is  to  be  dis- 
couraged as  its  greater  frequency  resolution  leads  to  greater  variability 
in  the  coefficients  of  the  Fourier  series  without  also  yielding  an  estimate 
of  the  confidence  to  be  placed  in  the  results.  Furthermore,  it  requires 
much  machine  time. 
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PROBLEM  153.QQ1C 


Direct  Spectral  Density  Computation 


Let  the  data  be  given  as  a set  of  readings 


Compute  the  mean  value 


and  the  deviation  from  the  mean 


We  compute  the  mean  lagged  product 


(16a) 


and  the  normalized  mean  lagged  products 


(16b) 


from  which  the  normalized  raw  spectral  density  estimates  are  computed 


and  the  normalized  refined  spectral  density  estimates  by 


(1)  Control  card  Information 


N number  of  data  points  (1150  max.) 
m maximum  number  of  delay  steps  (98  max.) 

At  interval  between  data  points  time 
C constants 
D constants 
E constants 
F constants 

(2)  Input  data 

Y.  data 
i 

i index  number 

(3)  Output 

(a)  Cr  mean  lagged  product 
(C j/Cq)  normalized  Cr's 
r 

rAt  delay  time 

(b)  V ' raw  spectral  density  estimates  using  normalized  Cr* s 

(c)  Ur'  refined  spectral  density  estimates  using  normalized  Cr's 
r 

Cr  If  C = (1/2  mAt),  then  Cr  = = r/2  mAt  . 

DUr  (See  last  line  below.) 

EUf  If  E = 2 , them  EUr  = lfr  « oee-eided  power  spectrum 

estimates. 


Prewhitened  Spectral  Density  Computation 


Let  the  data  be  given  as  a set  of  readings 


i = 0,  1,  2,  ...  N 


which  is  prewhitened  by 


Yi  = Yi  + bYi-l 


We  then  compute  the  mean  value 


-1  S b £ 1 


a = i Y_  ?. 

" iVl  * 

and  the  mean  lagged  product 


i = 1,  2,  ...  N 


T i i=H-r  _ 

' - IP — V Y • Y • , 
r N-r  / ii+r 

L i = l J 


a r = 0,  1,  . . . m . 


In  Part  II,  the  raw  spectral  density  estimate  is  formed  by 
m - 1 

V = C +2  Y C cos  ■aDr  + C (-l)r  r = 0,  1,  ...  m (17) 

r o q = i q mm 

and  the  refined  spectral  density  estimate,  corrected  for  prewhitening  by 


^ =T o— 1 r— 7 f -i  V,  ,+-?  v+4  V..1  r = 0,  1 

r 1+  b2  + 2b  cos-^  L4  r'l  2 r 4 r+lj 
• m 


, 2 , ...  m 


which  follows  from  Eq.  (24),  (26)  and  (27). 

In  the  computation,  b is  set  at  - 0.6  , so  the  refined  spectral 
density  estimate  is  computed  as 


jl.36  -1.2  cos  J 


* 

[K-.4vKn 


r = 1,  2,  ...  m - 1 
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where 


Pq  t n-i 


1.36  - 1.20  cos 


1 1 1 
— V +—  v 


X 2 *o  2 *1 


-x  V +-i  v 

1.36  -1.20  cos  7T- gjj[|  2 m_1  2 m 


However,  the  value  of  b can  be  altered  by  inserting  both  b and  2b 
into  the  appropriate  location  in  the  program. 

(1)  Control  card  information 

N number  of  data  points 
-b.  constant 


(2)  Input  data 


Y.  i = 1,  2,  3,  ...  N 


i.  index  number 


(3)  Output 


(a)  first  card 


Y mean  value 


(b)  second  card,  discard 

(c)  remaining  cards 


q index  number,  (q  = k-1) 


Y -7  fluctuation  from  mean 

q 

Y + bY  . 

q q-i 


PROGRAM  153,002  Part  I 

(1)  Control  card  information 

m (fixed  point) 
m,  (floating  point) 

N 

(2)  Input  data  — output  of  153.00 

(3)  Output 

(a)  First  card 

m (fixed  point) 
m (floating  point) 

N 

(b)  Succeeding  cards 
C 

r 

r,  (fixed  point) 
r 

C^.  (fixed  point) 
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PROGRAM  153.002  Part  II 


( 1 ) Control  card 

2 

b constant,  into  location 
2b  constant,  into  location 
<2>  Input  data,  output  of  Part  I 
(3)  Output  data 
(a)  V 


(b)  CO' 
o r 


) 


constants 


Cr,  if  C - 2m&t  < t^*en  Cr  = f 


r 2mAt 


DU  see  FU 
r r 

EU  if  E = 2 , then  EU  =W 
r r r 

FUr  if  D and  F are  set  equal  to  two  times 
the  confidence  limits  given  in  Table  I,  then  the 
range  within  which  the  true  value  of  the  power 
density  lies  is  given  directly. 
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